{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Project: Transient Navier-Stokes Equations\n",
    "\n",
    "The purpose of this project is to implement numerical methods for solving\n",
    "time-dependent Navier-Stokes equations in two dimensions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test problem: Flow past a cylinder\n",
    "\n",
    "The domain is \n",
    "\n",
    "$$ (0,2.2)\\times (0,0.41) - \\{(x,y) | (x-0.2)^2+(y-0.2)^2 <= 0.05^2\\}. $$\n",
    "\n",
    "The center of the cylinder is slightly off the center of the channel\n",
    "vertically which eventually leads to asymmetry in the flow. \n",
    "\n",
    "We are solving the Navier-Stokes equation\n",
    "\n",
    "$$ \\boldsymbol u_t -\\nu \\Delta \\boldsymbol u + (\\boldsymbol u\\cdot \\nabla) \\boldsymbol u + \\nabla p = f $$\n",
    "\n",
    "$$ \\nabla \\cdot \\boldsymbol u = 0 $$\n",
    "\n",
    "The time-dependent inflow boundary condition on the left is\n",
    "\n",
    "$$ \\boldsymbol u(t,0,y) = 0.41^{-2}\\sin(\\pi t/8)(6y(0.41-y), 0) $$\n",
    "\n",
    "The outflow boundary condition on the right $$ \\{x= 2.2\\} $$ is\n",
    "\n",
    "$$ \\nu \\partial_n \\boldsymbol u - p \\boldsymbol n = 0 $$\n",
    "\n",
    "On the other part of the boundary, no-slip boundary condition $$ \\boldsymbol u = 0 $$\n",
    "is imposed.\n",
    "\n",
    "The initial condition is zero $$ \\boldsymbol u(0; x,y) = \\boldsymbol 0. $$\n",
    "\n",
    "The body force is zero too $$ \\boldsymbol f = \\boldsymbol 0. $$\n",
    "\n",
    "We choose $$ \\nu = 10^{-3}. $$ Since the maximum velocity is one and the\n",
    "diameter of the cylinder is 0.1, the Reynolds number of the flow is 100."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: Mesh\n",
    "\n",
    "Download the mesh [flowpastcylindermesh.mat](http://math.uci.edu/~chenlong/226/flowpastcylindermesh.mat) and load it in Matlab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA\nB3RJTUUH4QgVDDMQpr0C2wAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAyMS1BdWctMjAxNyAyMDo1MToxNrb7SqgAACAA\nSURBVHic7b1BjBtHlq4bfu8CqhWZmzEEg/nABtyAAXpBvF7Ixc3Qi3enZtUFqMDqVaO6+wroZRVq\n1k/lu6dVsxtAjTYbs3FxJKC8Y6MX5mxI6wHX4ADihRY2zEEShtDeZHElzabf4s88FYzMDEYmZStU\n+j8YDTWLDEYmI+OPc+LEOe/87W9/U4QQQsjr5v943R0ghBBClKIgEUII8QQKEiGEEC+gIBFCCPEC\nChIhhBAvoCARQgjxAgoSIYQQL6AgEUII8QIKEiGEEC+gIBFCCPECChIhhBAvoCARQgjxAgoSIYQQ\nL6AgEUII8QIKEiGEEC+gIBFCCPECChIhhBAvoCARQgjxAgoSIYQQL6AgEUII8QIKEiGEEC+gIBFC\nCPECChIhhBAvoCARQgjxAgoSIYQQL6AgEUII8QIKEiGEEC+gIBFCCPECChIhhBAvoCARQgjxAgoS\nIYQQL6AgEUII8QIKEiGEEC+gIBFCCPECChIhhBAvoCARQgjxAgoSIYQQL6AgEUII8QIKEiGEEC+g\nIBFCCPECChIhhBAvoCARQgjxAgoSIYQQL6AgEUII8QIKEiGEEC+gIBFCCPECChIhhBAvoCARQgjx\nAgoSIYQQL6AgEUII8QIKEiGEEC+gIBFCCPECChIhhBAvoCARQgjxAgoSIYQQL6AgEUII8QIKEiGE\nEC+gIBFCCPECChIhhBAvoCARQgjxAgoSIYQQL6AgEUII8QIKEiGEEC+gIBFCCPECChIhhBAvoCAR\nQgjxAgoSIYQQL6AgEUII8QIKEiGEEC+gIBFCCPECChIhhBAvoCARQgjxAgoSIYQQL6AgEUII8QIK\nEiGEEC+gIBFCCPECChIhhBAvoCARQgjxAgoSIYQQL6AgEUII8QIKEiGEEC+gIBFCCPECChIhhBAv\noCARQgjxAgoSIYQQL6AgEUII8QIKEiGEEC+gIBFCCPECChIhhBAvoCARQgjxAgoSIYQQL6AgEUII\n8QIKEiGEEC+gIBFCCPECChIhhBAvoCARQgjxAgoSIYQQL6AgEUII8QIKEiGEEC+gIBFCCPECChIh\nhBAvoCARQgjxAgoSIYQQL6AgEUII8QIKEiGEEC+gIBFCCPECChIhhBAvoCARQgjxAgoSIYQQL6Ag\nEUII8QIKEiGEEC+gIBFCCPECChIhhBAvoCARQgjxAgoSIYQQL6AgEUII8QIKEiGEEC+gIBFCCPEC\nChIhhBAvoCARQgjxAgoSIYQQL6AgEUII8QIKEiGEEC+gIBFCCPECChIhhBAvoCARQgjxAgoSIYQQ\nL6AgEUII8QIKEiGEEC+gIBFCCPECChIhhBAvoCARQgjxAgoSIYQQL6AgEUII8QIKEiGEEC+gIBFC\nCPECChIhhBAvoCARQgjxAgoSIYQQL6AgEUII8QIKEiGEEC+gIBFCCPECChIhhBAvoCARQgjxAgoS\nIYQQL6AgEUII8QIKEiGEEC+gIBFCCPECChIhhBAvoCARQgjxAgoSIYQQL6AgEUII8QIKEiGEEC+g\nIBFCCPECChIhhBAvoCARQgjxAgoSIYQQL6AgEUII8QIKEiGEEC+gIBFCCPGC//a6O/Aj8sknn7zu\nLhBCiEfcv3//dXfBxo0VpE8++eTs7KzdawVhbeObx/1pENbavZb9PUqpogYXk+ViGhX9dTacx9Gq\ne7praby5GzY7DUsH8BUb+2k0W+r923+j8e1KKctVFxFHK9zMCh+cDefVvhRsecdU+ls7DjwLuPkb\nR4V7U1tel2B/EMryqm4XeLVX+gr79ko6tv2dX0yWym9NurGCpJRq91oYT81OaH/nuD8NGvUgrBcN\nl8UkCsI5Jsrj83vZN1xGo2BZU0p1TzvZv86G8yCsNXfD3J6M+xOMsNzPSgcw1SqlitrJfmlzN1RK\nWa7Lzvnw4f753uxi7viNxrfjyZkN50ePDks9QotJFEcry90o+tSgP9w/3xv3JxU6rJQa3B22ey3c\ntwofV+nPhJu2f75XoQUQR6txf4praff2tpwQ8TtWvi06eBC6p51xf5L7IJRtbTacHz3qXZ6Mtm9N\nbTdiDfAToG9Hp+UGcJbBZKuRCWbD5M7v9/aqNTJWk2pf/ZNxk/eQgrC2f7437k8Xk8jyNkydGC5F\n7xn3p5he273W5fEo9z3tXmsxjbLfhQmue9rBAsdgMYnG/enRo8N4eVX07YtJdHkywgS9/2Dv8iS/\nA5k+T7qnu7iuOFq5fEQHs3O71+qe7jp+o9HnZifsnna6p53BwYX9JzCIl6uyz/+4P7k8GR096rV7\nraBRL9nZpAWl1P75Hn7HCi0opcb96f6DPch/qUs2uDwedU93cf8tw9KtS5PmbtjutYpGYLmOnYxw\ngc3dcMuOSWvNTmh5rNzZcsSafTse7Z/vYQxv2TcsJdu91tGjw8uTUYWHEWB5ffTocNyfbn/z/eQm\nC5LSNMny+2HqtIjNYhLFy6vuaQdWVBytsq0tplG792H3tJN9EmYXczzAuZKDVZjlEhaTaHAw3H+w\np5SCtdfcDTc+IbPhPGjU5brKDl+5ZKVUtfkC7iaVPoc/6iM0uDtcTJbHT+5h2djsNMoqymISLSbL\no8c9pRSMpGrdaHYa6MM2c2Lqdewopdq9D+NoVVnbsNyBrQaZ3OZX0Bf43dPOYrLcUnSlNVxmtdsO\nth+xRmsqvWNYiW553+BDDsIa1mfVuoRVGua02XB+IzXphguSSn+/xWRZ9PtpU+eHuSbOuD/V3UcY\nDcbDE0crTP1Bo67/Cc8JhMH4k0pnB/wVjRhfHUerwcHw6FGv2QnjKNGz7mln4yQlz4CqNKlh6Sr/\nt+jOWMANwb9/vEcojlaQAWhJ+nV1+ModgQEqHrZmJwwa9bJTLd4v4wSNVJhh42iFtTn+bxDWttE2\nY7mDX6HaCl3Xtu07tphEi2mkt1bZlAeyaANbKtzlyUjfhtxGfWVpmHasVU0s4+UK05RSKghrR48O\n1XbLCz+5+YKkrJoURyuZOoOwlnWI4cnBWqnZacTRFUbDuD+RASqLF6WU4fqD2xf/7p7uzi7WtGrc\nn8pMGjTqIjnSt/M7D6FG+L/4RxDW2octy1xgPANl5w64PnQnNS7ZvQV4QfVX5BHa3jMjLCbR+Z2H\n3dNdY7cpaNQs/k8DSP7+g7V9mmanUWouQyNGJIXdCVwE1Gi9M2E1/xguwfgdKzugsqb8NoaIseJR\n21kPg7vD7umucaWVFS6704PWqnnb4H3VX4FYlv1BjckhCGvt3odKqfM7D8t2yWfeCkFSxZpkRHNl\nH359VdjcDbH0FunCADUWL/pTqk/NQVjX7YzsEx4vr4c7VsqiRuhq0EjmKZhiRWNaN4/063KZO3TX\nh459/yzbSHbTFY9QENbO7zysvBAWZNMo74tK7CEZNxnAIizVCPYb1rtRc7znAn7QbARK97QzG87L\nWqi6pSVgoJY1HcSUN16vYDqrdWfdemtVtqYM21Sopr6IZXhVrY37k2xQlcxFZX/QbDvd006713ol\nD5QnvC2CpAo0KTt16g+/bh4ppYKwLktvLA+xoIPZJC2Ii+zyeKRPLpjTMdNln3CYX/g3ZhNj0afW\np9oi94thHunX5eK4yy5djetymcvEC2ogj9Dg4GKbR8jYNMp+i3KLKdB3fYwW3L12RSqi4OqZRo5X\nKpF12T/hvpUKScD4ycaGVDAd7B0ru1FvOOsMKjjHsrapgEFYSuGyFqpQdjOpSNtUJZOr6J0SN3Qz\nNOktEiSVjgPdXpbYaP098vAbz6Ex02FuvTwexdFKlwq41Mb9qe6vA9gwz33C9Z2PXDXSd2VU8ZIt\nax7pvbK73bLOOqMFx7nM6KpB97Rz9OhwcHBR9GwbAm+0nN00yumqQ6Advr0ouBy/4MZGDL+r2Y0y\nESWw1YouHILnaNkUGQ0VeqWsc3SV1opXPEqbqR1byzrrsq25G5cSDlf0hlJ6iSfRft/cTa54eZW7\nyFMSv3c8ugFbSm+XIKl0HldKjfsTqFF2xGBEXh6PssczjZlOjIagUcOhCmzdx9EVzsHhW/BILCZR\nENaDRn1wcJGdesQdN7g7bB+aqrCY5NgcWfdLkXkk12Vx9BU569Y66eC4yG4g5bZz9OjQEmlS1MPc\nTaMsGwPtkhMwxapWFBVpsDFI0tGpVeQT03HflDL25HN75RjnsnGOVmVM5yJnnY77TK2HI1pay419\nzaVoMae35mjZYIDZ+9Y97bibXPZF3sbQrTeFt06QVLqZAbeYSoeOaAn+w4nO7HNYNNNdnoxw4Eke\n8iCsYQzhxXF/enkyOr/zED6ccX96eTyCUGFwwx+I5X/2e3OP5ojJIq9sfKIscVb2pauwMXo4dwMp\nS9lHyLJpVNCNwkA7CWSwd29jmNzg7jAIaxun143eNruZpTflMlPDUbyxV46GSHZPvqi1jaaz3Vmn\n47LtL25th9actqYcT606biYVOTkNYHJt1HK7GknHboAmvUWChEXcuD8Z3B0ODi4wCBbTSIREXzDi\n558N5+d3Hg7uDuVhk7gGPDCDgwvJ7LJ/vof/uqed5m4YR6vmbtjcDfHi0ePe8ZN7MIzgPWt2QsxE\n53ce4j98Re6qSg/k09EnKbt5JO/PfZzszjoDuxukaAMptzN4aAd3h/Z32jeNstg7AAt1Y1NGVKQB\nDEqXSQedsZgjG80sYaO9hUnfJcmFi7zZHWJGaxunaccVj3JzteW6tYvY6Gqz7JNl2biZhLHhkh7F\nVcvdnil5oF5hLOtPzE0WJGjG5fFocHd49l5/cHCB+aV92Np/sHf2/SmeyWanoWsJ/gvCOrx5+w/2\n2octRGCf33k4G87j5RWkSCl1/OQeHn4jYhXWlRFsnQYH48BjC+khjh73zr4/leQ68fLq8niU++QU\nBY/JJLXRPAJZR5+Ls269JzY3iMtqzmiq2WkURQo5bhplmq0XOdwGd4fd047bJFsvmv1dbCytHdvu\nHcxxx4k192TCemtTx16pTfJWdlTYTWcXZ52OfYy5OOsyrdlOPuRGJFqwK1ypX2GbePei1l5VLOtP\nz00WpNlwvpgsMRbPvj89fnLv6HEPkQhysgcDK/sUYXJv7oaSWAjag3NLs+FcpAi/uuEWgPvYOGIp\n+WDQiP51eFwx3QRh7fJkZMiSZYmUDuih2rQJIe83HH3uS1ehyA3isoGUxcgwJEEi7ptGBkUHjWGK\nOfZQj4o0KLU8V+nuXdYzI7nvHNtR6ZHb3Hm/lLapTfJWYVQUmTXuzjodjLHsYr8oot2O5dSUyz6Z\ngWUzKXv8ayMbT8s6usGFVxLL+lq4yYIkIyz3t0x3bvJjwRH1j7BdvBikmXox2envD9LUqLIXJRl5\nJVhLj+lq7ob6MTdYAHjAgrCOIDTI0vmdh9ocfW126O7Hs/f6GMpxtDp7rw8fIzaoipweurumlLNO\nJ3eRWPbJEbIZhspuGhlkDxpjpihlaeWmEbKH5xXRPd3NqkiFSV8VbASWcjoJRQfUqqUBLTJrql2m\nSuPmjTFmj/qzUBR84ehaMCgb5urSt+K409K6IrGs2yRk+un5P8/Ozl53H34U/v3f//35zrJ72nk+\n/+vlyeiDvZ/v1G/pb3g+/2scrdqHrZ36rWYn/Orh18jxo5Qa96e3W+9+sPf+Tv3W8/kPz0bffLD3\nvlIKrt4P9t7fqe3s1G+hWaXUbDj/6N4vduq3Ptj7ORZNH/zj+7db7yqldmo7z0bfxMur2XD++7/8\nGl8dL6/iZXLMDTF1H937hVLq2ejb262/C8I6uvTB3s+DsD66/+Xo/hgfHPen408no/vjZ3/+BhLV\n7IQf/Y//u9kJX1y9VO+o3//l1+3eh7dbf7dT30F+ttm/zUf3x7Ph/Nno28U0erF6+Xz+w4url83O\n/zX+dPJi9XIxjX71x/0Kd3infguqic6D0f0vu6cd41a7N9jshKP7X8ZLzBrv/P4vvy51ylVHv8lK\nqcUkGt3/Un4CR4Kw/mz0bRDWpBuLSXR58uejx4dlrxFNqXcUBoZS6vJ4FIR1/e45slO/tVPfGfen\n7cPrRf3nv/1i/3yvwu263XoX2eiNayx7r6S15/Mfns//Kne+8mUqpfBMff67L+T5nQ3nz+c/7P3P\nj6u11uyExmywTfdw62TeUOmip9pvillI/yGE8aeTCo8V7t64P5UeLqZR89b73W63bPd+Mm5y+Ql1\n7VGtDw4ukPBY/qTHrSUG+PFosRs1O+FiGiHPjdI2h+BdOX5ybzGJxpNkKTo4uJBA8DhazYZPlVKL\naRQvr2Q/PF5ejftm+Phisox7q8vjUfuwVeQrCMJaU4XN3XAWzdF++7AVNHJ8Mud3Hu4/SJLFyeJR\nbzaOVnF0FS9XYm+l53OncEsGjVoQ1suuOsUNIgvzUhtIRvcQfyiLQWynNTth7iW7sJgs1WnSPlIC\nVmgEaYSkAwhAqHCNSqn98z2MQ5V6sY6fVKy50NwNZxdzcZBWcBMJssUlnXEPssgFdkPyKG13mUo7\n54QThJfHoy1bQ/dgKOOhPvv+tHKD+ryhlNqme5iFBgcX2XIt1R4rvU1V3qZ/LdxwQQJwRs+GT8/v\nPLwOH1g/fSmaFAznQaOuaxV8GijcorQN8+5pp937ED/22Xv9IE2umhRhSrd8xIOBrFOYQTDbWvYh\nIG/j/hRbWXIh2Xdqib3rRVujCO3TW0CYBl6HixIHp3Dt7jLQ7n2I7a5mJ3TfQMLXxdGVFC5TafDS\n0aOe+PShoHG0WhwMg7SglGPHJB5SFeQHckR+YpWspiuqo9LGErYfqnmxpCmRt+2n6XavNbuYj/sT\nlNtRVbVN+oYIguMn97a8zLR7H+LU52KyrLwa0FrTL7a0k9NAVxEkCtqmexLgoP+aW+4D4czfbPj0\njQgHfysESWnnYcVUyuaGwttwWlZ3DWP2bO6GMKWDNK9MENZhEsXRSjbeB5OkyJscv4VEYX0Hhzgm\n3KBRx4SbzcggUnT85F4Q1s6HD/cfJIWdsn2eXcwlub3Md/a7IWlbYSfJ0kk3pEQG7BKlP5CWDaSs\nAgWNerPTaB+2svsB6FVamKCl0mcSd8+5Y8m6oSg/kCOBlkZoy8W+Sg3uyvt2RsckUYglI4Aj+BGD\nsD7uTytcY3bSDBr1s/f6jqnzjN2+LEFYw17sYpocBFzLjdLQV5ZrOVNyW8PFSjiS/as3Ij/E9sND\nJWVFr/SnuFr15EwPP0xmlbN/2LKHPypviyABsWngGoqjK7gUZJZMjoxMo/3D6zl9MYnUVMXLK5gU\n7V4rXl6N+9PFNOqe7iJOaTFZirvv7PHp4O4QfjB4+QZ3h5is22FLf/YgTrL215uFFOFt0M79sG4s\nnVTGVyNZrSzTE6otiMUgpf9UgSGF59+iBOJU0U/AXDsw1xWoe7rb7JR2ByW9CltGxxaTaHYxz7Xt\n8BEsLLZ0ViAyZTGNch1Z1z5GSUW4XKn1STZJwpv+bxyt4uWV5eiuXV30gQE/Z7PT2H79GzTqiPac\nDZ/qArP276V5UbkdDhrJagBDOvsG86vLq6ke6VC2h3Cz49ErOgNXKmG8fIXY0NV84AAWYbWA1Wzf\n9Mdwy9Z+At4uQVJaqGvqpJq3e632YUuCr2bDOU5EyiuLSdQ93YW0YKtDKQXZ6J52FpMoaNSx5kJp\nTqUtwWQKg9lxeTKCNy8Ia4tpdPS4J5bT7CKJm9KlSGmB1LIQ0w0gMY/k6hCNXWQkJbV/HiSpqZud\nEClkLNnA2mFSpkxplorKSBRenA2f6uoO701ZBQrCzYlN0TFdn3TbLj5Jwg5n0Tw3vE3lLeqLEH/m\nuD9V/Wky1eZNdpjpdH8v/pHUDWnU4mXSlO7UzfnGpa1vWal7hcTLK6Ua+tZ6WfsD4LyXvqm5DYO7\nw6NHPWhbheVFJiLxCo70oFEviohzDw+B0xtTClaWuT5wpZSjt1m8Dni/e+QqLtNYCOqPof9eu7dI\nkMQVhukg3Xpp6OP7uty4msjUj21P/MD6HvJishyrCaYV6NzgYIg1V+pVqF2ejBKvxcU8Xl7Bnz7u\nT4/P78nWSxDW4uFVvLzC8JUdV6APR5h38kruVrbFSDLUKG2ztZi6jnixVNS6RM2GTyVKAqd9XVqz\nUHaeNWy7OFph48cy6bt7zOQxTl2j6Uxdfp7FxBovV9jDKPvx3NbG/WkQ1rdfTZ/1+/vne7OLeYWA\nbwM9A7oxniu3Jk6CavHo+v+dDZ8GjTqioivbMTp4rOLlCuWhC9dJUitgN1kOFkmUrJvtPkA9Ggj/\nG6SnU6q5Il47N1yQxGIVHRL7A5bQYhrpkQ5K5tx0H1Wlz5XYQIODIcIf9s/3ZsOnECG9TFb3dFep\n5Fg4hikcRyiZg1RjCG2aDedBWMdf4dnLdkn3gxnDdHYx1wN/5T25RhJi8LLju91rjftTVTXOCLtu\niVk5nMN3F0clTvjntNnYaoKAGiXHKk9G3dNtZ0OZvNR2u/36bpYeI1e5NSy9ZZreZmLF2EskZNM8\naAcnotACogG3UV+9tSBNDHi0xa+AfF2YB/SIu8roZ7aMnzXXB67S3dAiNzg0UkJYDU84FCgbDbTf\nMytyvYncZEGaDecIa252GrmRnUFY7+IwUBrpIMaHBN3hhLl4kCQN2uxirpcvgvvu7L0+nFeLyRL+\nimYnTCaO3eQjqesvGZfZTAQSfAFNMmydZlpCEOmoi8rwGEZSkRqpNEFO5bhSuOBxRQjlgJafDx9m\no1d/ApBeXaxAxCNUflAl1l8pteVUaByn1UPAq3UMjaj1wOjKrcXLKyh3u9dCKa/KremlmMT7VNnq\nMqpy6DFy1bqnB+4b0fMVgF7K9CIXWzTyc3dDleEGT319CsfpopXuCYdfJ92OfeMVyOAmC5KCudBp\noNZvEXr0NpxX+ghQ6RoZb8aJFoyqwcGwuRt2T3cXk+ViskR8XbvXQhYf2C7ndx6KzgWNGswmBEdI\nqojcLjV3Q31Tyvgr9jwtBW90IwnxskVrXrx5MY3giHNHDBGjG0FYO3rcQ7q/7F9dsGSis/dHbE15\nEdnZqx5jinRbAZNXNXnThQ24h0TmtjY4GOqtGb7cshjB2clBoqpXqta3eYI0fUMFqys3Bj3V8g+r\nuUz1kMtAi54v2xQwfsEKi4OsG1xpvj4M6VflCfefG546COdbBwcXRRlLQRDWjh4dIvzp8niUniPb\n3T/fg7dXRgO2B+XYLGKlYJgvptFsOL88GQVpHb9k7ZME8kWDgyFeOXrUO35yD7nAcS4VOSDi5ZW2\nP1HHnxbTaHB3KBnBz+88HBxcwKaZXcwvj0dGGQsgtUovj0eLydI+lO1prXOZDeeDg4v9B3trZ43X\nO4C7JNmPflSQ9S6bgxX3sEKDcSZ9auBQ4dCxKZCbF8cFBEmuhY1tSrpqIZslKNiUitTet2yYgGMB\nCIOiqhyBWwGI3AZVJiYicCvqkUtuEjz3WlMWEFSCHSmV+NUn/ockbM/Nt5BgAM2GTxFfkD3/oQc7\nqHSAYtQik4JSanBwAaul3WshyBs2U7IBM70elHjFcLXF0WqhouMn987vPEQisnF/2uw0MOZwbA1O\nNiiTBGLh48ZeOvoDP2G8TILXoWf4UjnYe37nYXM3dCi0U84iwdkXY8GbXa7KzYephyTE7t/iDtx0\nuUdfg7BuDyPMRc5pGZ+Cv6ish6foWG61HZGiM0ySdLWUSYpJP+vQzubgcMGS0ltP3+DYmiVhRIUf\nAsuCXCtNP9zt2BrAL2u8KFbXNvtwWL/un+8h5ggzEo72V/M6vCncZAtJwJN/9Oiw2QlRfwhrjTi6\nGvcniEc4fnIP+5xYySJ5uxydg3cF2yRyNAfuY3F5SawE1GgxjSBaScRdI9l5xjYSChLn7t5jxwu7\nMhKNg2A8/JdUYNtNzo12TztSb+n4yb2jR4dwLienXpZXG02EQDv7aSeOVtgSyxU58UMaL2IzyVKz\n3KWdov7Iz1cUrVSUsduCJbODe9lWgD22opluY6kkA3teV0v1xeIGLZO+ax1YYM9cXtbq2pgwwqWM\nkI6lQnw129dSLGobq0tpJwWT8x7pWQJMYovJ8g0tLeHCWyFIAKPk6HEPRgxWTHG0Ovv+VBbv+O1R\nDQERDbOL+dl7/USBlldy4lLaFNNEP/YoAXL6RpG8U6ViNjgYLqaRBCjD7yfTDR5vQyrE7ZBbuzYV\nrbqcT+qedmYX841+s2anYa9qqlI3XbViEBWepY1vwzaV4TbM0u61LEdQs+SWkBdkf86lKdxzi4Va\nappebCq7LhsYLq2pTRnwArfaccLGDHjN4gIQBi4ldEs57jbmfELSL/dbt7FYlGPp+ixIcq8vsIKw\nJufS8KOgtMSN9OC9RYKk0pJ9MGvEssndXgoaNXixcAAW2yFwFMBwwbEh/QS+8UVq/UQ9TCL8G2Ga\nSFKH5wrbCUePDiVrkR4dLlIRJyX+dlXx7gj2jdqHLagsXHbYVbbIUrv3od1rhyB4WJlF78lWfFj7\na1jDasDlWbI792CoOdaQtdTZy1JUQl4HdUY2TtPidbG/rVlQAMKgaCPKwH0DI3aoKuQ+6Ttmd3W0\nuhyzuzrmJXKsxlTKvtyYo6/arp6MaqMp/bGS5Z1SSpw9N4a3QpCgQwgHUErBIYZwcMOPBxtlcHeY\nHCDdDTFMU29bmkdnGrkMXJErMaRU6sFTCNLrNJCTEQdyMdQk39118aRUKvR8rNgdyX4j4sjF2MdQ\nRnwHZClfgIu9duKmM1JIVMMx2MGibUXxC4VNpSuPje90zDPkaIhkQw+KkAgUe99c8sO6m1x6cLaF\ntrUOrNHaxi91sbrcs7u6tOYo5KqM+joWi5ITGhsbVOlTVjSqcz3hN9KDd5PrIS3+61vUE7o8GTU7\n4Uf/4xd7//NjGUbP539VSn2w9/7t1rvtw9YHez//6uHXz/78zWK6VO8oHGTBLLZT21lMl/Fy9WL1\n8nbrXdSteLF6+WL10rEzeOeL1cud+q2d+q0grC+mS5XkE1q+WL1MKirVdpRSfZvvrgAAIABJREFU\ncbQa3R8joZFSaqd+C3V0ns9/eD7/QRZ6O/VbaYDfdRmbf/nv/6pPW4tptFPfSSoz1W/dbr37wd7P\n4+XV+NPp8/kPt1vv6hVWjAJCepv753vtQ1voPECZho05V6Ta0+e/+yLbDXtT4/5k/On0V3/8pUt/\nhCCszYb/O3uI2Gj5+fwHZ5GrGzWEDAZ3hx/sve/YSdQ3mg3nKLuV21oQ1hyr7KDwkl6kJ8tiEi2m\nS8eAhWYnHN3/8nbrXcvFtg9bpS7WKOakd+zy5M/H/59rOMBO/daL1UupWJYFZaIcoxVQvczysyql\n4mj1+W+/cCwWdbv1blGJI52NT9mz0be5twuPklJqdP/LF6uXGy/T/3pIN9lCgt0ju0TGr6WXJ1Bp\n8hvZScJnF9No3J8mezla8IKxHmk2m2dnZ2dnZ19++eV333333XffffbZZ0dHR9kupYl25qnXrtPs\nNPAidkSkwiPOzKbnjXYR4W3MIPo2kpEyFRjGvtIWVhJlINeS9QFm3dl2dGe3y5tLBTtsjF+wfteG\nMEIUM3Q/52G4Ug0qlJRN7O88k1E/A+vIRu9TbnB2EXKQKPevG3dTslisrgqlmCxuwLKlNIK0MLnl\nPaUiD13CJXBWxG7+WkbvDfPg3WRBEnB8R87r4EWZpMb9ydl7fYljPnrUi5dX++d7cthIpc4cCa7T\nGz87O/vuu+/u379///79brfbbDabzebR0dFnn3323XffZQ1QOZgNR+J1tu/dEDFy2N/qnnb2H+wF\nYW12MUfsA+RQn7NEQjCmszN1rltPrQ9i0QN8Vir4De4OcU9+pFhtZXU7GNrmGL9g+aKgUS/abEi0\nvOSpw6IAuY2hB0U9zHW1OW5E5bXWKZqbLMHZRVgOElWuwr6YLI27V60UU5HjrvIPYdnSyz14ZMce\nLnF5PJpdzO1rLJc0rzfGg3eTXXbxuz/s1Ha6p7sf3ftFczdEVe+v/vD1V3/4+tno22ejb57Pf0BF\n5N//5deoZT66P0aF2cuTEWbzZifcqe2od1QcrbI+ui+//DLXEgJBEHS73W63+6c//en6xbAG3x2c\nfr//y69nw/neJx9/9Yevb7feHd0f733y8YvVyyCsf7D3frMTtg9bi8kSpa/hzUOf8XE8Ic9G3+b6\nEF5cvSwy9pXmOkOV9xerlzu1HVwjHAhlKzE/n//VxW+Q7Ub7sPVi9VJ3OzwbfaPSZOSf//aLOFpt\nU85cKRWEta8efp29FZJwtmzjuRXc42j1L//9X3/1x19W6Gquq21wMPzVH395+8N3y7a2U9vJdSri\nETh6fFj0wSJyvU/j/mSntlOtYrd6R43ufymfLeusM1ozHHfb/BDZmu7Cv/w//1qhTjx8nnrddJW6\n/nbqt3712S/tH5d5aeMXyaP01cOvc322dNm9TtKl0zSOrhByun++Jyd1JBRbpScrZQUNNxdWSRLI\nEKR5EoUvv/zS5aftdrufffaZ/F9ZvGBRDDdUvFwhbSvymuiuNjjog0Yd/T/7/hSBGEg4ppQqcscr\ntxOvsraC32M2nOP84E+cJssIdgjSbBel4hcs5N4K2fGudrE4i6obXtuUplUZV9s2dfyKTK5xf1qt\nhGu2QRxdqHxIU473Sce2qZtuhBfayzHbKfKzWQ4ebWzQCJdA2q1mp+Fo+7qfzFNK4WCiejM9eDfZ\nQlr817cwMkb3x/p6B2EFXz38uvtPnefzv+4/2IPxNBvOX6xePvvzNzu1HaXeUe+oF6uXH937xbPR\nN7Bp9DFxdnZmsY0M2u32f/7nf85mM3kFrSG04cXq5U7tVhytELkAk0UpFYT1z3/7RfuwhaOdRnjC\n7da7z+c/wHJ6Pv9h/OkEZhPeAyRcIrdXsPm++sP/Wkyjz3/7xfP5DyqNv3j252+yrW0kXq6ez38o\n2l7eiBHs8OLqv2b/Ni8bv2Bp/NnoW30YSNzzNtLb7IRiJCUh41v0Vl/pwwe1TYHtIKw/n//wfP5X\nuUC02f2nihIShHWMKPzEn//2i+4/mVuzpbjdenf86QQWiVKqgqUlYPB8/rsvPrr3i+1v3e3Wu4a1\nuphEs3+b/+qP+xW7V9t5NvoGT+tiEn3+uy/2HzgFCoFno29vt/7O3TKTR+mrh1+PP52Icea/hXTD\nUwcpLW93V3NP675gHMFr91qDu0NUukwTdTeUUrOLuW5LCffv3y/Vjfv37w8GA/m/aK17uit1MbB0\nlWoUSqnFZKkFedeMUrODg6FegzxOq/xJTCriJvSDuiqvkCtuAhanlycjlCxDdXa0lmQX3rpGjp1Y\nq3Gun+5Cl+RUo2OJs1wQhiB5erZZRAuy5YA7vH1OF2Ragsm+fT1s5H1HHtJYq+OwTfdwbAA7fNsk\nyZYeorTd0aNe3mmEElXyMJgHd4fx8mr7W2ekcK22VaZ3Ly3aeTUbziskwncPF9K/FGmOL49HRuE3\nb7n5gqTyNEmSUOn55xfT6Oz708UkQqY4vGIccQXutpHQbDa73e54PL5+ZTeU6hiSQRwiAb0sWrzr\nqdtQE0ydqiCtoGooUxDWkvNV1kKukkD68niE4O/uaUedXlcKvxyO0GdIXf59bpRwLGSlEcrXPmxB\np+PlFZx4eokzpZRePwYWj2O5aH1b3sj6XA30p9kJx/1JPFzB4abWK7oWfcr+OgYGliZi1SHRVKmC\ncuIsOnrcgzm4ZYiKNrGu9LIs9qrtKo0T0/3VxjtzQxYdUywatzQIa+d3HkrGHeMGuoeMSsVbx4NH\n0hmMWLkbiOZFlhYpxeTSlNDsNOyDykL3tBP3VsiDp5Tq/u4fqrXz0/BWCJLSHqSjR4frx0uT1xeT\nSGZzfESKzmXjof/+7/++Qh8MQQLx8mo2TJ5VycqKP437E/QKJ3Lw8Ev9oeS6MtnwdGXCxBE3VlkF\nMpBTnMbQl9ZUasRA6lS6B1DqKU2LvtikEeHdiTRG17madBWMjUKcaTpatV6L05i75fDvxkrY6WS6\nNq3Ekh0qzf8khwTiNK27/l1FjRem6kl/yqS6dlhrdsI4SrJSyVer9Kh1+qnkGrO6pVIDDmcGsvcw\n/ceV2qQo2artar0si/TBuEy5KOmPcdWz4RzFlF9JdR+YR0Gjvn++h58PFyWxo9nfLhk2Wm9l2Eh1\nQT3/bHZsJD9QXsvJ0D1sYTTCl1i5LEtlkhMmuyEq4/jMTRYklBjXX0GRoThaNTsNY7tPTv+0ey0U\nI4caxdFqFpnhwtWcsIaMISmRSlxzyVpsMY3G/QnkBBVpMT3hlWyu36KwBWyQBI06ctnZR//g7lDO\nacFuyy0gq0sdPA/jyTQ+Gelmk3QmLihtqTYVV5YkmM1OA8djc3tiUSmV1v6QzlyX4wxrMEb1qGiL\n0ujzLEROZWwUWaxYbEd3MMvjJltak2lRaVoi065xLXG0UtPEbsjaKGJJqE2KIn9FFnncxi2vF0kX\n4+VqyzK1SjIHJoXJzRAkHePuWRRLKTXuR0FYG9wdSnCTWtcbGRgu51LlrIUq4+C1PJWOSDqY6k38\nJNxkQcol/1xOY+28jr7vrdIh+MpD+5NBli7qFco6LK+On9zDeakgLd+AhHsqze298Yow+LqnHcxr\n9kgbY+3sGJgnmzrXPr3Uv4RH13DBOS5+pSy3SndTXD6l9wr/LirHKeKBCciuNI5gSg1eRR1xKQaB\ndbTlnVn9sLSJdfGrOlWGVVr3UefyeDRW1Su3qvUMeNtXbsWuapCeObM05XL3MHLwK8TRCqFrW949\nMW1Rceb8zk9RVRkHG4JG/fjJPf+D7m6yIBn7eHFa5AZPlAThwFeWrEkPWwg1VkqN+1NUg1VaAjqw\nWCyazeaW3dPbxPwoNpMOamSIKBpGEmZhGeiYyDD4rt9QXMYbloQeVB2k8dburnaYTVI3HSnVnW/D\nWmeUlpLA3nNHgrQcJ4Q5yYrbqL+S1b1an1K3rCOutNDnV3LtSk5Zne/B7N6ySI+0iXgZ7MuWKnFk\nIDnp1dY13fWASexBbmmtYmnV7rWwtptdzOOoXEIKg/UlbzL8HN13lWsoz4ZPZ8N55YMNPz03+RyS\nTqyVXEti0taPLHRPO/HyCoMjkFzgk+ViGiFLwvZ9WCwW8m+0r9LEEHh4jp/cW0yWuitJ0rzuP9gL\nGvWjxz0cGNJPGEgeUklnYFZNzStUoYozFMjGQClw/h8n8Ct8PNYSmQtFPS+L3Jn2YUs6+UpWi/qU\numWpUKNEgiU1kSPXZ34bNThv270WFltbNZtWSAnSRDvVnAdGgvDAmiLBjmxqpk2VyO+eCx69JGV+\nox40amldzepjBtsB+is4e+fYbNmbPBvOK6faeo28FYIUrxcADbSDb1JeCGvSy+PR+Z2H+w/22r1W\ns9NIQ3RqxvJEz7zgjv4p2YfH8zw4GAaN+mz4tNlpLCZLFAPE0TljPGEW0GUpXl4h8R18MtnBZ2Tt\nkw4UZUGuIANwtaHie7y8qlKWO+9IaeUC5EK8XqgCyo2ubjm/qLwptXLl72yJBD3Le7W+STYppANW\nSkGTKheOQ7OytabSlX4pz+p1U5qWg2o13bPFk/BDVB452H+VlPnSphzcribAuWFyjs2KI8Tti1aD\nu8Nxf3L0qPdGhHrr3HxBivPKUUuCqcuTEdx0mNmRMiSJZRjOsaiRTW8hGyzngv4psZDg/VBKtQ9b\nEq6jlFpMI0vqNl2WxOVYXKwzx9635BTIFTALokboGLwc7h9XxUnMpAB5qdaEbKKHOFrJogSJv7bR\npNnF3MiR0XSuQWeQXRzIXkiFjkmGtOyfsPbappipccntXqvCJefWT0qDTcrZhbkpHrDzWqodEKcF\nXPR1BsJG0j3dFlzTFVouitBB4kp7lmF7sTEBdXayC9k3hRsuSLlqFKdnRVE3FrKklMKMgHgzESFM\nzUbQzmKx+M1vflOqJx9//LHRsTS/6hR1+eLoSnZNMa1fnoyu/XIZUYEbSuKnYaHnzl+yISGv2E/h\nlHJYYybStbDsXGBJghk4lzLKgp/VskgUEa2mSchynd2oqOC4K8pJgyzvZTuGrLj6/TQmsm00KTdS\nC5dc6jZmzSNQtqZ7USnY7Jh3IXe6UOs3EH421Kss1bjd89/shFuukMQZgOIG1Rp57dxkQZLhhUM8\nOLF89l4flfrkqIFKa5Y3O2FzNwwaySQIKwpzohGGp5Qaj8fudpLxZtmjwr/h8YCiJG633RADVGUS\nUuHJP3uvj71WvLl9mNTfQ7XybOZjhFAnnXEojuC4AIRr0dCSUpsBFs8hkLgSd+LierKLaaSf3ArC\nWrv3oUpltRRFadzKOu7wuxQ0VXovBNZqzr7geihXNU2SsinZxnOzd1vaUQVB0qXunr0UbNlNOJyI\nz6pR0fWWdfnqp8dykRVSrvsuCGtFKzPMCYODi+7p7vZZH18vN1mQ4H8bHAwHBxdy2PDoUQ8VYzEF\noLw37FwJdcO4xBvkfJLS4rNVaiTp2YCKGI/Hhnmk0hkf7g7oXxytsN8u7wnStKez4TzNozqRjcrj\nJ/cMHzfC2+Sgg14ZVqwWl5z8aUDHhsklV40ANgPsH5dG7NlIy8YXbczHamSjwU1GSFWpb8k1j0DT\nuVRobK0jjt/U8U7GacnRbGtxtMrm4EFSnFJTqiXqr1SAQ5F5BJCy1qVj9qUMNuEcl1YSAOIeXFpq\nS0nyBW9stsgCK/q4RAW/oT46g5ssSNgGbO6G+w/2oEB6mT4ZRnp5b1TkU2lxP5kL4EzDzCgjY7FY\nfPLJJ5988okePqeDN2TVCMBqGfensMZwhCV3WhTFEnvcMrhlTCNiQkqzK6WQyduxnLN9HrSokXLe\n/3Cpf1PK97LRTVe0SkWlRPcItI1ZrrunHReLQU4BF73BcUMOITAIFt34ZiDWobsm2U9WOgY4WMwj\nYWOBQeWQezsdhE/t/VGb1MheVMxxSwnp/Df2RFreP9+L05pkRd2QKNw3MXihiJuc7fvFz1f75//4\nYvVi/On0qz98jSJDUpLk2Z+/2andQt7infqtndpOHK2Q8Xqnfuuje79AMZUXVy+bncbo/lgKkOvD\nIo7j8Xj8xRdf/Md//Eccx3Ec48V//ud//tOf/vSb3/wm69Zr7oZp+fN3grAehPW9Tz5+Pv/hg398\n/6N7v4iXV5cnI6Xe2anfgjp+9Yf/hdRBSYHzTJmTZ6NvgrBupOVGut+P7v3iduvdxXQ5uv9lHK2e\njb5xLI7wfP7XxWRZVNXCpdp3ENb0ajdZFpNo/OnUqRT0Oyq3lJFOHK0+/+0XSqmjx4eWjJwWIWl2\nwherl5cnI6NuTVHPf/VHWxmbbL2fvM5szr2NIkmBtQY2CmD/6o+/LPplv3r49Qd772dbQNr75/O/\nLqabzxLNhvOd2i37r3C79a6eDjyX0f87xjliSzsbC5OjJNjGI19BWBt/OrUnEZ8N56P7Xx4/uVfU\nJXtRMSTVHhxc2CuBoZSUuwWDjP4vVi+keLzRjXF/glIAv/rjvnsKWv+zfd9kC0lp5sL+gz2sIgd3\nh/qSU7ZksKuEwGWccMTKUV+9os509lsWi8VgMPjNb37z8ccf/+xnP/vZz352dna20ZuXJNs+bMEy\nQ+NBWA8adUjRYrJEtAzmUPHGlHIu6dkHEEDoYm1YYo5nw7lLtW84Xizf5Wir2TsjXYJ9YJ+hNrpN\nYENvXPCO+1OX06/2UqHu9Uz1/b/cdlwKYBfNWbCTXPbSZxf5aZwM7AEOLuaRtFO0f2b3c+psHITI\npLflYWG47zZuKZUt6ydz17g/Hfcn4rvWt0hvjGEk3HBBAlibIE66fdiSnX9MZCrdkumeduLoCvkd\n4CBOq54vz74/lRiEbQ7JBmkiLPHRzS7miOhDqMXsYo6yF1BHI20dxiicSzJpxtEqm2JV4STK8Qih\nd9iswo705clI314q6Gf+w4MH2HHj1LI7XarWmd1rlx65ONzYmovbRE85kfsGbKQ7ZgGwuJ7cyxlY\nohZx2GjLc/gSkmCfTx2v2h7gYN89MtqBFz37p41+Th2LnI/7E8fB7JJMK9m7LVgsVj5cr0c6qDS2\nFlukjnfgzeKtECQBctLsNGKtSlAcrWV0xlyMlBsqlZDzOw9xSknOCZUdDdigxm4Q/h0vr5q7IWrX\nItri+Mm9o8e9IKzLtCgf1wc0Ti0UBZ6KzYetI4ngWEyjdu9DY3upaJYP0gRC+oul1EgVr09dwvwy\nTeVMK1BxGJGv8OHEniISwGf/ihzwjk0FmVKhoFQp2CI9hhHvspW9cYdjoyaVyjJXFODgbh4BmJjG\n7+6y77jeSL6cXx6PXAx95WzZ5C4WhY0hdkUYsxN8zjfPMBJuci47A2R2QoY6VDeQSnR6Rm0pr7D/\nYA8FDlD3ReL0ZJ9ZL5hkIMeqJYJ8NpzPormcuMbjiqo/WClLtuNxf4I1bxxdiWvCGNDNTrgf1nXF\nkqtDDyVbvv4GOeR0nbP1ZKTSA/zmJawnECqrRsCoiadSV1VZJ0lzN0RBJgEGbinjII6uHCcFLHgv\nj0ft5VrKbRyXzt5bC+1ey0gbKslBSjSSuY0IKnG/jRsvPEhz0+XmS80eAd7YGlZUeg/dzSPBSHBX\nYfAEeblWYcf8GBHSaZWHtQx1jiF2QpyWCsMKWF/LYjglKeEPW9vUq/STt0KQYDTg57yua5IaOlhx\nwDGNxVQSjd0JkXtUyrMODoZHjw6xFx2kxVhVGgolY07P7JAEW2OLSAsUbvf2cKxV6vcgUBh7SJLu\nGmZQczfM2vtBWNt/sJdmHkqsriJbIftIBGGtHSajHEPcSPKIBELoCU5xVXC146L0cOFqlTfF2MKZ\nZdgcZftTdFQ+F5mg9ZSa4/6k7JSq0lkVHto4TUpdqgXjNpaaT91zClg0aTGNSq9Fei19RVXWPJIu\n4UwbGqk2eIxcq3piEcc+lMrLgBOEeoGJjb7iuKCOs77eGkyGWM4iAA8bbJh2SlUm85wbLkhYRyul\n2r3W8bk5f+nG+Gz4NF5eyVAY96fxyQgVzVFKNUgrumIQYH6B3ohXIV6vpBJHKzGngrCGzSEYMUl+\n8UZSpVuleduw8Em7J3PimndRLzKEE7tY0VtWYUWPRJDm6oYBdD68TiUuJVhgHVbb+A3SZDBH6Uzq\n7qrKNoW7jRysP4HXwpigK5hH0o7Uu9KTUpdqAaEuzU5Ydj5VZdzLuZpUuSoE3JXj/iQtylVay1Va\nhR0J66oNHjlfXO3uVSANc3g6uDs8etzLLgfTRC1X9mKVOotpdPSoh0cy0FLsY0LAfGUv6PxGcJMF\nCSeK9s/37L8QJlyYF4O7w2StkU7EcVoWBSNGap8j7RAS+4ulsn++h1wvGGR4thFKAFchvhHvRCyf\nUqp72sNiR3a2sps3QVg7e68fpEWGZODiT3hoLf6rjYUMYBTClESFi3bvQ6RJ3bJymqzukRCslJBo\ngRtXSikc23KMXM8iaardCdKCwipdXlT4XpXWHoVlU3Vy3708GSEWv1QLcXRVKn17VpMwwkv3WCl1\nLSdLVcY8MiwSpNFSSrV7e8aaz4UgzbU67k/L3j1poewmULIKURMUUw+SDJlX+lJSf5DtpDkyzEPi\ngVQmO10r6Ix17ZbVuV4LN1mQgjSBd+L1Wl87pOns8rci4jRjqRRME0+dqJT82PAkIKpNYpTFs7SY\nRvFwpdJRhVw44/4EOYqw5MFJI732qNaTK+xU5a7NpdAc9jyaBXOuy0ETpZ3IQzCP+JcMgcxu81rG\nPYxOnP+FsMWZep3plZqVW6VlmKp48fKnXQnqrlE9i2iRG6eogK/8++y9/jb9mV3Msbxw3D+QVN/u\niCah1pGjq1N+VuM3Ven6zDh3bPeDZU2KIE2yVfTZ3MuEGGPwVJ6dkdi07Md1vz2qcLkrkIGsCSw3\nTcymWKuZ+cY59G6yIKUmbVLZenYxT1Ixpj8PTAF9+R+ENQQyYCJudkJYQio9pt7sNPQSrlL5WIrV\nI7kWHBSShSEI6ygoDr8cFi+Xx6Nmp5Hk9k5Pk4iDUbokCcRyz14gKbjatCmNjBWO9w3PHkRUpYbm\n2hsyc272OdGfXtFvTEmiMSpd4uFteuXWbAuSH0UhBjKNPXF/3tynpKxHRZIZXieBzgu1V3lqjWvB\nEYLK5t24Pxn3p7jSxSSaXczlV8CLr3Z/+9pdfDLCHKdSsVGymEiT0+vrs2zd93avhacDgdH6Vzh2\nRpZ6dldb7mSdlgqb4hEwJoFXe9OyHnU4Hsb9KYqZVWsWnnz0Ezm9NiY3KXLoxdFK/UO1XvxE3GRB\nAjLlYZbHXCPr3MVkeRmN9KGJR10vnqTSow8qPeo/PpjACbB/vnd+5yEW0c1OCK+aPAMKJ1Eu5u3D\ntUjuZifU94cQTIFNoGA9dk5pFbIxoxlmgb56LdIkFy8HtlWxTYIHCT2Pl9c5yN3JxPs+hblZeeMH\n9//yZBRHV81O2D3twEcRR1ez4Xw8uXagF00xG/elZWMZm3zGehbiZ0nluRF81giWcwencWcXyZEy\n6TNEIom+ORii53DjWGZbw0iVf2t/ui79Povm+Iexhmi6lX7HSi5eXuG3q3Dts4s5fOP2t+VbSGEN\nR8dw4B0GHy4NNw1bxc1Oo5oZkV27wAsiZpCEO1Wu/6uX/NBDjTZiOPSSoywUJN+YXcwRXXb06NB4\nnvGGdq+FFxHchSDOpCbF8gqLNZxvh92ThPEsV5hzZ8OneIDF14eoWXFezS7m2GqS0kdGxXERntlw\nLrmC9GCt5EIy+hQkeYifYicZL1qONGZ1SByDl8ej9mELBdTLzsKGcSOTqSodEKCUdvTE2Au5ft7S\neUHijvI8tFeGjZi7nlVK5RoxcbTCkWr9xrqTLHIzMYeOIAYym7DVWGyp9D6otDh9fJJExAzuDrOO\nUJUxaFRqAkJpxv3E22YPmbGzmETNTiMIW9U2IyWPrfgqy7aAfcfZcC5ljZRSsCGUNnLGkykUPbuy\nCdKSSMC+dskCLay8F4u8/tIUnPxleVM2k94uQYIvDjFaSYyA9jzDBMG0JS4RWTDKBg8cL/Br4d/w\n0mCqwpod74RLUCVbONe5O7Ecw+swpPBODDtsXyeCdHE9FoPM2Y7c5zMIa+3ehxLho/LO31h06Ppe\nYTtqN0Rah8qeDZz3csx+nQXhVbhkY17QETeFBMXqzhmsMLA1LSc88EFjPVvYjWmkHTH5sOzjjUVu\nkB5iKy9I0+s4T+tvgfug0iGNAa+UgpFdqtvJOgA5fw8uKmvSYhp1O53seSxHZDrOHsZy/LhSqtlB\neuL8mq26OMkKVV/ZqETgzaA4Fwcsfi+c63DMeGT0v62V/Cib/F7vxmKyrBbo+FPyFgkS4heKxhCC\n3PCD6QGvsIfkTBL0CSqCMPHBwTBIY7IlBwECPbFjJMtSPNLXgXxpBPNisoT1E0er7mknSWc3iRCu\nqrue270W9ioxrIvCxvQIH1SVFTNiNnyqGwRFQcwzrVJ1tYkAYHnbPe0ppbLnE13A8lZW8S5FM/Xd\nXfGh6wfCxMXq3g3Zf5IA7lJXAe+uUqrd+7Bswe/0BHeyQ+PusVGJh3m3ggTKEVSJ1DdOuboj+5f6\nSssduXUV1jTwh6PbLraFvkI1jvug/5ayJkWIf0Li10v9FtkQx1KHopI+JFvUhy65z18vb0XqoLi4\nYhuQNAR6jneMZkgUBuh1mp/HvaPHPXhRjp/c23+wl7jXG3UpUxSEdayMlEqymqJZeIFVms8N0dXi\nFRkcXMyGT9uHLeRqzK5okARIwsQti1bJjY/3D+4Oz+88RMJWvSJULmtu65JFPHXG/ak8ThXqn5bN\nE2MAh173tINrkb0xvQqJC3rZBewMl7obugyLS9bxs7Bi5ZdqH7bcyxXOtNJ/pTocr1dNjJerdqU6\n5SrdQdH2nJwKHQk4NYF/l711ytDy8rYFBgyeMpjXFby1cZqNLChfoB2LYH2sBm61yowO4Nd8I7x2\nN1+QNlZsM5LiiCBhIx2CgVcwvUo5VIwV3Y9kZLJCy5JBfHBwgcdJHIDvvnhHAAAQqElEQVRS2QgW\nFRLNSZjDbDhHxRodSNflyWgxMXeGMFFCycb9CTQYVh18hmffn6Ig4eabpqXOC0qWQNVvrNLD0krW\nP83mw27uhmWrxyp4Yi/mcBuqqgXLNZ9J6buBHSD5v93TXfc+GE6ejYnPhVhLiQ27yr3DetVEuXCM\n1Qp3T9+6Q7ZZ988aPoBSa5rUOjeje0qBg7TN3bDa5RvxRGXXdrmJE2VF68jGGpheccMFaWPFNkON\nJJYXxTehNyqdHJVSWHHr/h/NuqoheymMEiOXcFL79WIOkdAnGn2vPkiSfO8mZ56OR4O7w8vjkfw3\n7k9mw6dYKqIdmD6ooHF5MpIntn3YQkUG8TU53jTdXwewti01lajMZJos7pyfJd26qoxeSBA/UFxc\nHKEIw83SdC4Iq1JXsPFxx5KD2cziMho3fhY3X5yu7kKO07u585djoQod49YFZcrbwzDVx2GpNY2R\nZwhWToV68DKGy9YzVJl4orKrmdwlaallmcxjju9/7dxkQRr3p7Ph3FKYIDdhaLy8MtJRQzAk94k8\nVLJ5bmzzwMpRmWxjMrgX00g8CWp9yZNE5Q3nqYtgF2Ua5D/5CLra7ITd092jR4fwwiFfOLLzwfWP\nsIhS/hbdXye4xN3q6N4SARtmLh9HqnLjhyvrdUHBzTSyYyX1O0rPqpnaAd3TzsytshSS8xovYltu\n42dzCy+5mDvGaTZ3uwpXpI9bfd8ujeHcUPhnrcG8W+coKohx118JnGu6y1n4tY+XyViRTTIUONeO\nErJbnhLdsPGzRjiD1g3Xp6BCWv3Xzk0WJDidYT2c33l4fuchrA1JcIdNGt2QT2tPrHuKOo0k3Vw6\nvjEfBWFt3bqqpyeQJs3da3d57hnvxWR5fudhEsC2G8JsQn4EpaUNjZcrhC/LfxCbIKwnWVkRolrs\nHcYarZRlMMurVF1qbQuDMvskOK5w8fEtE45hQsn6aQO38j86ucd+HfcDcm9mei7H5kGScGfj9Y3m\nTpypX+doVyEkz37b4TF2qc4unTEGp+OtSw4PZZaSbYea7nLSwHjd3dmFcZ6zGig5eHKTXECSN/oP\nLRURXXyP7hUgveImR9lJIHX2nHlSDbZRRwEIcfXKPzCdyYtxtNIEJs1wk7fJCRfN0eNeshW0e/0A\nSAKIcV9J5QsEgCGcwUgbgaOguVFJ12c7jkeWMSfb6UGak23j6b+idZkqEyaE0Ljs67LCtbdQlNT5\n+rfYtD1rz6EZWLNaGEhyJuN1HPu3343iRS4iMJ9avr3IY7nRW6g7666/0cE4gHN7472VsbQxELwo\nqNLlPFbRiqSp5X0v+qwlutplKtfdvPIpceaXGzx5cbBBWiXL8uTC+5L7BolrsNyBuFJSeR+4yRaS\nEKQnKGFkyIoVwXLHT+6dfX969v3p0aND/Gn/wV77sCUusjRnTAOfgpcMbUrxWUxbug8wWC9TJrtZ\nMGggQt3TztHjHgJM9x/sGRuwlp0beEKSI1PFC0Z9kRWkOdns98oy0B3XtvZCAxtXuPaM4BsjvxFR\n2T5sbVzpOy5142V+TlX5fS2fzTqdBPt9sNxDu7mTTT0Fmp2G/Yez3PZcA3H/wd7GWu9FYykorgYr\nWE4I2Gu6F90BlUZIWr5U5alRlhKDp2DxtDG6wQiEMT9uNfXiaDU4uHiDAhl03gpB0pGQ/OyfMNPJ\nQQRxkcFji0gb7MdIWN3R415S7PwkiYtrdhpwCc6Gc4l6GPenetB5ENZ0rwu8dvozhldU8c6NBPjZ\nt3aMPVVIst3tlutium7BIUzIXjQIgfVFLWQjo8yPW1f6sEpzMzrHmULvjtOKpbqa3Y1Z5HQC9tCG\nDfewYBsp66wT9GDR3K9TBZsNRfVSEd1jP1OV3UDS+pNTDVbvj2U6ttR0V+nmZe6fNu6+uKhR2tTm\nwWMfOfboBntmentcw+XxqOzBBn94uwRJD8kPMnW3cEw1WykVSV/WXryYS8ZJCdcOtDyhi0mE/65P\nxa6HL+sPBmIc9GdMMjTnnr0wzrUUaUzuGtO+mWTx18nX2R8kOUdsacGyL72xAltQnKxBzo25H710\n1CRLEWvLFv3GU8BFoQ16Ms2iD+bOR7nOOoDRntvPypsNWLEhKi+XTefkCsPfczcghaCgpruStCA2\n87rQQrKoUby8ys1wj8FTtKqw1+WzxGpmw1wzX12orNuUOPGBt0iQ7JYsjqka54pUulY1bGRJJCNj\ncTGNoGSwovbP9/CfHMnUNWO9MOBcvHMyOvU4iOxza3hCoDHZpyJ3U9T+FLnsD1kciThftTEYochb\n5VK+r8hlF0er8zsPLZORpUG7JlmW+crqxtxYfglB8NkfYmNd2tyoOYurSj6VfdE4A5slaNgm8e5p\np9lpuK+H1vpTcE42G+2d99kcrx1M0s33vEDJHG0jndRFMS16muxXgdio7O11KRif+6PgZr5xgQw6\nb5EgXR6PihLG6Gfa9SlPgkd1G1k/PS7SJUHhxtBMop8bdV2T9CWenG9Fipek2fR0t8oTAGOKLLJa\ninKqGptbxn1wqVJR5Cd0TNWV663C/60Wogo1sjvNLb4auyZtjKGQNLj6iy6zapCmttNfzI1Xzn5Q\nrY80i7NOyJ3Etz81WWRzuyR5yj0na9l40780a2Q7HlzLWtgodVZtHrdrksW2Vlp0g/G6JRuyfDCr\nrG9oWJ3B2yJI2OheO6OXRmnrKYJUusmBf6PYhFq3kWUV1uw0pJIpjgQZBz8Tq2s3jJdXRafqtNau\np+nFNNI3PAwjKTtFZheb9vVpkOZp1V/c6K/TP57dO0HPHX0FhsWJdbpL5ses93wxiRy3cC2zQ5Em\nOSbfyyq0JWZXRw5ZCzLk7BjbSBZnnZA9lfxKTk0W3brFZOnwi5gDyb7xpn/QWNM4ZpmSZ1aYDeeL\nyXKbebxohWe3rZP+ZDZlHYdc7hO0/UHy185bIUgIAi76mQ2/s/gosLMqD3mslenDizI5Sk7udu/a\nuS/783omK3luMZ6MfHSyqaBbSGrdSCoar1hsXmct2jQhZrOTlUr7mN07uTwp3EzOYgRx2F3/FhCi\nUjZTai4VzifpnzVv5qZFrnxQn1il4P3GD+rbSBuddel31eXEgnI+NelyDDPXD+wyHavMQHI/f6bv\nwBUdPMpiLGhyj8Zn2Wgl43iWEXbocj4h696wB3QIxoW8WfmBLNx8QSo64wZ0Zx2QJ1Bfq4qNvB5I\nndTVltlHf3plf14/3iRTHsaTkY9OPy9pDGUxkjaE0krhQYcJEc4WEQZHf53+dTIjuDiadHRzMPsT\nWL9XM1VfnRqljZua5C7S+k6e45wCUExLZfKo2pFtJBdnHcAYhn3wyt07YiWII1rlHd7K/aA+kNwH\nob6R5mIgCvIpRzVyxAiFL3EhWnSD/cjE+tddPwhvXH4gCzdckOSYau5f4+Was04nd60qJcmvW4hW\n+uwjC15jgg7SJFoy5cXLVby8MvZgxYORc7Y/NZIsC085luRo9UvUOKYnx7W58XWIJHR0NK19vNOQ\ndHwVyrRAjYrSt2dxTKxpaJLjMh8fbB+2sNR1vP/pB+vx8gpjJptsyfJ1CjehzFws9/zHcO8Eadrf\nOFoVHd7K79VuKHfAfRBKAGopX7E4Hl6tGgEJhc9dU1qQ6IbFZENQhg6+BWvQNys/kIWbnKlhMVnG\ny3lRERdEGWQ9RTBoZpH5QRwtzGYLnQ3X3ok8Q3G00ksNZdOnIqpH5VhCHSnrZwDT3uIHCNIj9EGj\nvjFKRz6CI4ql7Bvts7uXJ0nN3FJiptLKQJhP3b8avw5OepUtz+M+08lRfBevi4AadIO7Q3vQdvbr\n2r0W1jpF5anyP5h6cUtN/ZfDUVn3jnuSbNmbLDWc9HOypWQSxmW8vHL/FH7NcX+y5b5REfgtqt0B\njGpnV0ENdlXlcAw/udGCNI0sZx0MZ7oORq1RzArvX0zW6k6ilJH+TrjpjRMG8fJq3J82tU3IIKyp\nqQryEtrLnJvb7dyPGKDYRKm9kLQuRqVilOW/bpvPBmENdefKfmPZb0G6wgqfKltKDl7fCgXoYIqV\nSsoXJ+VIypUwqNAxVX444T64dwzvD8JaqU+ptG6hez66ssNAlb8Dci2lauiVfRAWk6X6B/fmXwPv\n/O1vf3vdffix+OSTT153FwghxCPu37//urtg4yYLEiGEkDeIGx7UQAgh5E2BgkQIIcQLKEiEEEK8\ngIJECCHECyhIhBBCvICCRAghxAsoSIQQQryAgkQIIcQLKEiEEEK8gIJECCHECyhIhBBCvICCRAgh\nxAsoSIQQQryAgkQIIcQLKEiEEEK8gIJECCHECyhIhBBCvICCRAghxAsoSIQQQryAgkQIIcQLKEiE\nEEK8gIJECCHECyhIhBBCvICCRAghxAsoSIQQQryAgkQIIcQLKEiEEEK8gIJECCHECyhIhBBCvICC\nRAghxAsoSIQQQryAgkQIIcQLKEiEEEK8gIJECCHECyhIhBBCvICCRAghxAsoSIQQQryAgkQIIcQL\nKEiEEEK8gIJECCHECyhIhBBCvICCRAghxAsoSIQQQryAgkQIIcQLKEiEEEK8gIJECCHECyhIhBBC\nvICCRAghxAsoSIQQQryAgkQIIcQLKEiEEEK8gIJECCHECyhIhBBCvICCRAghxAsoSIQQQryAgkQI\nIcQLKEiEEEK8gIJECCHECyhIhBBCvICCRAghxAsoSIQQQryAgkQIIcQLKEiEEEK8gIJECCHECyhI\nhBBCvICCRAghxAsoSIQQQryAgkQIIcQLKEiEEEK8gIJECCHECyhIhBBCvICCRAghxAsoSIQQQryA\ngkQIIcQLKEiEEEK8gIJECCHECyhIhBBCvICCRAghxAsoSIQQQryAgkQIIcQLKEiEEEK8gIJECCHE\nCyhIhBBCvICCRAghxAsoSIQQQryAgkQIIcQLKEiEEEK8gIJECCHECyhIhBBCvICCRAghxAsoSIQQ\nQryAgkQIIcQLKEiEEEK8gIJECCHECyhIhBBCvICCRAghxAsoSIQQQryAgkQIIcQLKEiEEEK8gIJE\nCCHECyhIhBBCvICCRAghxAsoSIQQQryAgkQIIcQLKEiEEEK8gIJECCHECyhIhBBCvICCRAghxAso\nSIQQQryAgkQIIcQLKEiEEEK8gIJECCHECyhIhBBCvICCRAghxAsoSIQQQryAgkQIIcQLKEiEEEK8\ngIJECCHECyhIhBBCvICCRAghxAsoSIQQQryAgkQIIcQLKEiEEEK8gIJECCHECyhIhBBCvICCRAgh\nxAsoSIQQQryAgkQIIcQLKEiEEEK8gIJECCHECyhIhBBCvICCRAghxAsoSIQQQryAgkQIIcQLKEiE\nEEK8gIJECCHECyhIhBBCvICCRAghxAsoSIQQQryAgkQIIcQLKEiEEEK8gIJECCHECyhIhBBCvICC\nRAghxAsoSIQQQryAgkQIIcQLKEiEEEK8gIJECCHECyhIhBBCvICCRAghxAsoSIQQQryAgkQIIcQL\nKEiEEEK8gIJECCHECyhIhBBCvICCRAghxAsoSIQQQryAgkQIIcQLKEiEEEK8gIJECCHECyhIhBBC\nvICCRAghxAsoSIQQQryAgkQIIcQLKEiEEEK8gIJECCHECyhIhBBCvICCRAghxAsoSIQQQryAgkQI\nIcQLKEiEEEK8gIJECCHECyhIhBBCvICCRAghxAsoSIQQQrzg/wck3/aymHuhwAAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "load flowpastcylindermesh\n",
    "showmesh(node,elem);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Uniform refine this coarse grid, and project the circular\n",
    "part back to the circle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: FEM for convection-diffusion-reaction equations\n",
    "\n",
    "Choose either non-conforming P1 - P0 or iso P2-P0 elements.\n",
    "\n",
    "Besides the matrix for Laplacian and divergence operator, you need to\n",
    "compute one more matrix for the convection term $$ ((\\boldsymbol w\\cdot \\nabla)\\boldsymbol u, \\boldsymbol \\phi) $$\n",
    "and the mass matrix for the reaction term $$ ( \\boldsymbol u, \\boldsymbol \\phi). $$\n",
    "\n",
    "The mass matrix for non-conforming P1 is diagonal. For P1 element, use\n",
    "mass lumping to compute a diagonal mass matrix.\n",
    "\n",
    "For convection term, first write out component-wise weakformulation and\n",
    "compute the element-wise entry and finally assemble the matrix. Note that\n",
    "the derivative of a linear basis is constant which can be factor out the\n",
    "integral and one-point quadrature is good enough.\n",
    "\n",
    "Test your code for a scalar convection-diffusion-reaction equation\n",
    "\n",
    "$$ - \\Delta u + w_1 \\partial _x u + w_2 \\partial _y u + \\alpha u = f $$\n",
    "\n",
    "by choosing a smooth function and moderate convection coefficient."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Time Discretization and Nonlinear Iteration\n",
    "\n",
    "* Implement Crack-Nicolson method for time discretization with `dt = 0.0025`  which leads in each discrete time step to a non-linear system of equations.\n",
    "\n",
    "* Use Picard (fixed point) iteration to solve the nonlinear problem. The fixed point iteration was stopped if the Euclidean norm of the residual vector was less than `1e-8`.\n",
    "\n",
    "* Discretization of the linear systems in space by either non-conforming P1 - P0 or isoP2-P0 elements. Solve the algebraic system using direct solver (backslash) in Matlab."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4: GMRES method for solving the linearied problem\n",
    "\n",
    "Suppose the linear saddle point system system is $\\begin{pmatrix}F & B'\\\\ B & 0 \\end{pmatrix}$. Use the block-triangular preconditioner $\\begin{pmatrix}F & B'\\\\ 0 & -BF^{-1}B'\\end{pmatrix}^{-1}$ and `gmres` to solve the linear system to the tolerence `1e-6`. \n",
    "\n",
    "The inverse of the Schur complement is approximated by the\n",
    "least-square commutator\n",
    "\n",
    "$$ (BF^{-1}B')^{-1} \\approx A_p^{-1}F_pA_p^{-1} $$\n",
    "\n",
    "where \n",
    "\n",
    "$$ F_p = BM_u^{-1}FM_h^{-1}B', \\quad A_p = BM_u^{-1}B'. $$\n",
    "\n",
    "The Poisson solver $A_p^{-1}$ can be replaced by the direct solver or\n",
    "one V-cycle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5: Evaluation of the Computational Results"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Matlab",
   "language": "matlab",
   "name": "matlab"
  },
  "language_info": {
   "codemirror_mode": "octave",
   "file_extension": ".m",
   "help_links": [
    {
     "text": "MetaKernel Magics",
     "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
    }
   ],
   "mimetype": "text/x-matlab",
   "name": "matlab",
   "version": "0.14.3"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "138px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
